NICE Online
2023-03-01
1 Setup
1.1 Libraries used
Load libraries necessary to run the script.
library(tidyverse) # general
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(ggdist) # plot distributions
library(ggside) # plot distributions on the plot margins
## Registered S3 method overwritten by 'ggside':
## method from
## +.gg ggplot2
library(easystats) # estimate and process statistical models
## # Attaching packages: easystats 0.6.0 (red = needs update)
## ✔ bayestestR 0.13.0 ✔ correlation 0.8.3
## ✔ datawizard 0.6.5 ✖ effectsize 0.8.2
## ✖ insight 0.18.8 ✖ modelbased 0.8.5
## ✖ performance 0.10.1 ✖ parameters 0.20.1
## ✖ report 0.5.5 ✔ see 0.7.4
##
## Restart the R-Session and update packages in red with `easystats::easystats_update()`.
library(patchwork) # combine multiple plots
library(modelsummary) # for tables
##
## Attaching package: 'modelsummary'
##
## The following object is masked from 'package:parameters':
##
## supported_models
##
## The following object is masked from 'package:insight':
##
## supported_models1.2 Reusable functions
1.2.1 Plot distributions
plot_distributions <- function(df, col) {
p1 <- df[c("Participant", col)] |>
group_by(Participant) |>
mean_qi(na.rm = TRUE) |>
rename("x" = all_of(col)) |>
mutate(Participant = fct_reorder(Participant, x)) |>
ggplot(aes(x = x, y = Participant)) +
ggdist::geom_pointinterval(aes(color = Participant, xmin = .lower, xmax = .upper), interval_size_range = c(0.03, 0.03), fatten_point = 1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
scale_color_material_d(guide = "none") +
theme_modern() +
theme(axis.text.y = element_blank(), axis.title.x = element_blank()) +
labs(y = "Participants") +
ggside::geom_xsidedensity(aes(y = stat(density)), fill = "#9C27B0", color = NA) +
ggside::geom_xsidedensity(data = df, aes_string(x = col, y = "stat(density)"), fill = "black", color = NA, alpha = 0.2) +
ggside::theme_ggside_void()
p2 <- df[c("Video", col)] |>
group_by(Video) |>
mean_qi(na.rm = TRUE) |>
rename("x" = all_of(col)) |>
mutate(Video = fct_reorder(Video, x)) |>
ggplot(aes(x = x, y = Video)) +
ggdist::geom_pointinterval(aes(color = Video, xmin = .lower, xmax = .upper), , interval_size_range = c(0.03, 0.03), fatten_point = 1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
scale_color_viridis_d(guide = "none") +
theme_modern() +
theme(axis.text.y = element_blank(), axis.title.x = element_blank()) +
labs(y = "Videos") +
ggside::geom_xsidedensity(aes(y = stat(density)), fill = col_green, color = NA) +
ggside::geom_xsidedensity(data = df, aes_string(x = col, y = "stat(density)"), fill = "black", color = NA, alpha = 0.2) +
ggside::theme_ggside_void()
wrap_elements(
p1 + p2 +
plot_annotation(
title = col,
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold"))
)
)
}2 Analysis
2.1 Load data
folder <- rstudioapi::selectDirectory(
caption = "Select Directory",
label = "Select"
)
df <- read.csv(paste0(folder, '/9_analysis/Online Study/data/data_pruned_for_analysis_2023-03-01.csv')) %>%
mutate(
Typology = fct_relevel(video_primary_category, "urban", "green"),
Valence = datawizard::rescale(Valence, to = c(1, 7), range = c(1, 9)),
Arousal = datawizard::rescale(Arousal, to = c(1, 7), range = c(1, 9))
) %>%
rename("Participant" = Anon_ID )
length(unique(df$Participant))
## [1] 377
video_summary <- df %>%
group_by(Video) %>%
tally() %>%
arrange(n) %>%
summarise(m = mean(n), s = sd(n), r = paste( range(n), collapse = " - "))
Each video was seen 23.16 times on average (SD = 1.876, range = 19 - 28).
2.2 Summary data per video (stimulus)
dfenv <- df %>%
group_by(video_name_df_keep, Typology) %>%
summarise_at(.vars = c( "Restoration", "WillingnessToWalk","Presence", "Width", "Structure", "Beauty","Interest", "Scenic", "Familiarity", "Valence","Arousal", "Crowdedness" ), .funs = function(x) mean(x, na.rm = T)) %>%
ungroup() %>%
mutate(Typology = recode_factor(Typology, "urban" = "Urban", "green" = "Green" )) %>%
rename(`Willingness to walk` = WillingnessToWalk)2.3 Participant
Here we only provide an overview to confirm we have loaded the correct amount of participants – as a sanity check. Participant characteristics are analysed and reported by a different file (Table 1).
dfsub <- df |>
group_by(Participant) |>
select(Participant, Age, Sex, Language, Crowd_Preference_Mean, starts_with("ipip"), SSA_Mean, SPS_Total_Mean, SIAS_Total_Mean, NSS, Upbringing_Environment, Current_Environment) |>
slice(1)
nrow(dfsub) # should be 377
## [1] 377A total of 377 participants (Mean age = 30.8, SD = 10.0, range: [18, 67]; Sex: 49.3% females, 50.7% males, 0.0% other) participants were included for analysis.
(dfsub |>
ggplot(aes(x = Age)) +
geom_density(fill = "orange") +
geom_vline(xintercept = mean(dfsub$Age), color = "red", linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme_modern()) |
(dfsub |>
ggplot(aes(x = Sex, fill = Sex)) +
geom_bar(stat = "count") +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_see_d(guide = "none") +
theme_modern())
2.4 Descriptives of rerponses
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
df %>%
select(Restoration:Arousal) %>%
datasummary(All(.) ~ (`Unique (#)` = N) + (`Missing (%)` = PercentMissing) + Mean + SD, data =.,
output = "latex_tabular", title = "Descriptive statistics of responses to the stimuli",
label = "response_descriptives") %>%
save_kable(file = "output/response_descriptives.tex")
## Warning: To compile a LaTeX document with this table, the following commands must be placed in the document preamble:
##
## \usepackage{booktabs}
## \usepackage{siunitx}
## \newcolumntype{d}{S[
## input-open-uncertainty=,
## input-close-uncertainty=,
## parse-numbers = false,
## table-align-text-pre=false,
## table-align-text-post=false
## ]}
##
## To disable `siunitx` and prevent `modelsummary` from wrapping numeric entries in `\num{}`, call:
##
## options("modelsummary_format_numeric_latex" = "plain")
## This warning appears once per session.3 Correlation plot
library(see)
library(correlation)
p_corrplot_u <- dfenv %>%
filter(Typology == "Urban") %>%
select(WLW = `Willingness to walk`, Crowdedness, Valence, Arousal, Restoration, Presence, Beauty, Interest, Structure, Width) %>%
correlation(., bayesian = T) %>%
summary(redundant = FALSE) %>%
plot() + scale_x_discrete(position = "top") + labs(title = "Urban spaces") +
theme_modern()
## Warning in genhypergeo_series_pos(U = c((n - 1)/2, (n - 1)/2), L = ((n + :
## Series not converged.
## Warning in genhypergeo_series_pos(U = c((n - 1)/2, (n - 1)/2), L = ((n + :
## Series not converged.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
p_corrplot_g <- dfenv %>%
filter(Typology == "Green") %>%
select(WLW = `Willingness to walk`, Crowdedness, Valence, Arousal, Restoration, Presence, Beauty, Interest, Structure, Width) %>%
correlation(., bayesian = T) %>%
summary(redundant = FALSE) %>%
plot() + scale_x_discrete(position = "top") + labs(title = "Green spaces") +
theme_modern()
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
p_corrplot_g
p_corr <- p_corrplot_u / p_corrplot_g + plot_annotation(title = "Correlation Matrix")
ggsave(p_corr, filename = "output/plots_corr.jpg", width = 13, height = 7, dpi = 300)
3.1 Video Ratings
3.1.1 Variable Distribution
3.2 Determinants of Perceived Crowdedness
3.2.1 Models
rm(preds, params, mods, vars)
## Warning in rm(preds, params, mods, vars): object 'preds' not found
## Warning in rm(preds, params, mods, vars): object 'params' not found
## Warning in rm(preds, params, mods, vars): object 'mods' not found
## Warning in rm(preds, params, mods, vars): object 'vars' not found
preds <- data.frame()
params <- data.frame()
mods <- list()
outcome_var <- "Crowdedness"
vars <- c("Beauty", "Scenic","Width", "Structure") # replaced familiarity with Beauty
i=1
for (i in 1:length(vars)) {
var <- vars[i]
print(var)
f <- paste0(outcome_var, " ~ log1p(video_mean_pedcounts) * ", var, " * Typology + Familiarity + (", var, "|Participant) + (1|Video)")
model <- glmmTMB::glmmTMB(as.formula(f), data = df)
mods[[i]] <- model
# Parameters
param <- parameters::parameters(model, effects = "fixed")
param$Variable <- var
params <- rbind(params, param)
# Predicted
pred <- estimate_relation(model, at = c("video_mean_pedcounts", var, "Typology"), length = c(50, 7))
names(pred)[names(pred) == var] <- "Value"
pred$Variable <- var
pred$Label <- paste0("Main effect: ", format_p(param$p[3]), "\nInteraction: ", format_p(param$p[4]))
preds <- rbind(preds, pred)
}
## [1] "Beauty"
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.4.1
## Current Matrix version is 1.5.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## [1] "Scenic"
## [1] "Width"
## [1] "Structure"3.2.2 plots
rm(p_crowd)
## Warning in rm(p_crowd): object 'p_crowd' not found
p_crowd <- preds |>
filter(Predicted < 7) |>
mutate(Value = as.factor(Value)) |>
ggplot(aes(x = video_mean_pedcounts, y = Predicted)) +
stat_density_2d(data = df, aes(y = Crowdedness, fill = ..density..), geom = "raster", contour = FALSE) +
geom_line(aes(color = Typology, size = Value, group = interaction(Typology, Value))) +
scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
scale_size_manual(values = seq(from = .05, to = 2, length.out= 7), breaks= c(1:7)) +
# scale_colour_viridis_b(option = "A") +
# scale_x_log10(expand = c(0, 0), breaks = seq(2, 20, by = 2)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), breaks = c(0, 1, 3, 5, 7)) +
labs(x = "Number of Pedestrians", y = "Perceived Crowdedness") +
# coord_cartesian(ylim = c(1, 7)) +
labs(alpha = "Rating", title = "Predictors of Crowdedness") +
theme_modern(axis.title.space = 10) +
theme(
strip.placement = "outside",
strip.text.x = element_text(hjust = 0, size = 12),
axis.title.y = element_text(face = "bold"),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
ggnewscale::new_scale_fill() +
scale_color_manual(values = c("urban" = col_urban, "green" = col_green )) + # col_green
scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
ggside::geom_xsidedensity(data = df, aes(y = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
ggside::geom_ysidedensity(data = df |>
select(all_of(c("Typology", vars))) |>
pivot_longer(-Typology, values_to = "Predicted", names_to = "Variable"), aes(x = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
ggside::theme_ggside_void() +
facet_wrap(~Variable, strip.position = "top", ncol = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
p_crowd
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
3.2.3 tables
library(modelsummary)
options("modelsummary_format_numeric_latex" = "plain")
names(mods) <- vars
modelsummary(mods, escape = FALSE,
stars = TRUE, statistic = "[{conf.low} - {conf.high}]",
coef_omit = "SD|Cor",
coef_rename = c(
"log1p(mean_pedcounts)" = "Ped. Counts",
"Typologygreen" = "Green",
"StructureTypologygreen" = "Structure x Green",
"BeautyTypologygreen" = "Beauty x Green",
"WidthTypologygreen" = "Width x Green",
"ScenicTypologygreen" = "Scenic x Green",
"log1p(mean_pedcounts)Typologygreen"= "Ped. Counts x Green",
"log1p(mean_pedcounts)Beauty" = "Ped. Counts x Beauty",
"log1p(mean_pedcounts)Scenic" = "Ped. Counts x Scenic",
"log1p(mean_pedcounts)Width" = "Ped. Counts x Width",
"log1p(mean_pedcounts)Structure" = "Ped. Counts x Structure",
"log1p(mean_pedcounts)ScenicTypologygreen" = "Ped. Counts x Scenic x Green",
"log1p(mean_pedcounts)WidthTypologygreen" = "Ped. Counts x Width x Green",
"log1p(mean_pedcounts)StructureTypologygreen" = "Ped. Counts x Structure x Green",
"log1p(mean_pedcounts)BeautyTypologygreen" = "Ped. Counts x Beauty x Green"
),
title = "Mixed effects models of Perceived Crowdedness \\label{tab:models_perceived_crowds}",
notes = "Ped. Counts : Log-transformed, average pedestrian count per stimulus",
output = "output/table_models_perceived_crowds.tex"
)3.3 Determinants of Beauty
rm(preds, mods, vars)
preds <- data.frame()
dat <- data.frame()
mods <- list()
vars <- c("Crowdedness", "Structure", "Scenic", "Width") # "Interest",
for (i in 1:length(vars)) {
var <- vars[i]
print(var)
f <- paste0("Beauty ~ poly(", var, ", 2) * Typology + Familiarity + (", var, "|Participant)")
model <- glmmTMB::glmmTMB(as.formula(f), data = na.omit(df[c(var, "Beauty", "Typology", "Familiarity", "Participant")]))
mods[[i]] <- model
# Individual-level correlations
# dfsub <- modelbased::estimate_grouplevel(model) |>
# modelbased::reshape_grouplevel() |>
# select(c(1, 3)) |>
# group_by(Participant) |>
# slice(1) |>
# data_rename(pattern = paste0("Participant_Coefficient_", var), paste0("Beauty_", var), regex = TRUE) |>
# merge(dfsub, by = "Participant")
# Parameters
param <- parameters::parameters(model, effects = "fixed")
# Predicted
pred <- estimate_relation(model, at = c(var, "Typology"), length = c(30))
names(pred)[names(pred) == var] <- "Value"
pred$Variable <- var
pred$Label <- paste0(
"Objective Crowding: ", format_p(param$p[2]),
"\n", var, ": ", format_p(param$p[3]),
"\nInteraction: ", format_p(param$p[4])
)
# Data for density
dat <- rbind(dat, data.frame(Value = df[[var]], Predicted = df$Beauty, Variable = var))
preds <- rbind(preds, pred)
}
## [1] "Crowdedness"
## [1] "Structure"
## [1] "Scenic"
## [1] "Width"3.3.1 plots
p_beauty <- preds |>
ggplot(aes(x = Value, y = Predicted)) +
stat_density_2d(data = dat, aes(fill = ..density..), geom = "raster", contour = FALSE) +
geom_line(aes(color = Typology), size = 1) +
scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(y = "Beauty", title = "Predictors of Beauty") +
theme_modern(axis.title.space = 10) +
theme(
axis.title.x = element_blank(),
strip.placement = "outside",
strip.text.x = element_text(hjust = 0.5, size = 12, face = "plain"),
axis.title.y = element_text(face = "bold"),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
ggnewscale::new_scale_fill() +
scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
ggside::geom_ysidedensity(data = mutate(df, Predicted = Beauty), aes(x = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
ggside::geom_xsidedensity(data = df |>
select(all_of(c("Typology", vars))) |>
pivot_longer(-Typology, values_to = "Value", names_to = "Variable"), aes(y = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
ggside::theme_ggside_void() +
facet_wrap(~Variable, strip.position = "top", scales = "free_x")
p_beauty
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
3.3.2 Tables
names(mods) <- vars
modelsummary(mods, stars = TRUE)| Crowdedness | Structure | Scenic | Width | |
|---|---|---|---|---|
| (Intercept) | 3.861*** | 4.006*** | 4.469*** | 4.071*** |
| (0.049) | (0.045) | (0.024) | (0.045) | |
| poly(Crowdedness, 2)1 | −0.972 | |||
| (2.226) | ||||
| poly(Crowdedness, 2)2 | −17.489*** | |||
| (1.568) | ||||
| Typologygreen | 1.328*** | 1.087*** | 0.243*** | 1.178*** |
| (0.042) | (0.040) | (0.028) | (0.036) | |
| Familiarity | 0.133*** | 0.113*** | 0.033*** | 0.099*** |
| (0.008) | (0.008) | (0.005) | (0.008) | |
| poly(Crowdedness, 2)1 × Typologygreen | −18.689*** | |||
| (4.482) | ||||
| poly(Crowdedness, 2)2 × Typologygreen | 9.603* | |||
| (4.064) | ||||
| poly(Structure, 2)1 | 55.409*** | |||
| (2.156) | ||||
| poly(Structure, 2)2 | 4.805** | |||
| (1.542) | ||||
| poly(Structure, 2)1 × Typologygreen | 0.514 | |||
| (4.849) | ||||
| poly(Structure, 2)2 × Typologygreen | −3.159 | |||
| (4.236) | ||||
| poly(Scenic, 2)1 | 137.104*** | |||
| (1.270) | ||||
| poly(Scenic, 2)2 | −2.104* | |||
| (0.963) | ||||
| poly(Scenic, 2)1 × Typologygreen | 1.391 | |||
| (3.208) | ||||
| poly(Scenic, 2)2 × Typologygreen | −4.448+ | |||
| (2.626) | ||||
| poly(Width, 2)1 | 59.713*** | |||
| (2.295) | ||||
| poly(Width, 2)2 | 5.092*** | |||
| (1.534) | ||||
| poly(Width, 2)1 × Typologygreen | −17.449*** | |||
| (3.576) | ||||
| poly(Width, 2)2 × Typologygreen | 1.021 | |||
| (3.480) | ||||
| SD (Intercept Participant) | 0.777 | 0.998 | 0.536 | 1.041 |
| SD (Crowdedness Participant) | 0.155 | |||
| Cor (Intercept~Crowdedness Participant) | −0.614 | |||
| SD (Structure Participant) | 0.166 | |||
| Cor (Intercept~Structure Participant) | −0.848 | |||
| SD (Scenic Participant) | 0.091 | |||
| Cor (Intercept~Scenic Participant) | −0.917 | |||
| SD (Width Participant) | 0.178 | |||
| Cor (Intercept~Width Participant) | −0.857 | |||
| SD (Observations) | 1.290 | 1.216 | 0.728 | 1.201 |
| Num.Obs. | 11107 | 11102 | 11104 | 11105 |
| R2 Marg. | 0.136 | 0.195 | 0.688 | 0.193 |
| R2 Cond. | 0.366 | 0.519 | 0.797 | 0.539 |
| AIC | 38293.3 | 36891.5 | 25248.7 | 36703.2 |
| BIC | 38373.8 | 36972.0 | 25329.1 | 36783.6 |
| ICC | 0.3 | 0.4 | 0.4 | 0.4 |
| RMSE | 1.26 | 1.19 | 0.71 | 1.17 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
modelsummary(mods, escape = FALSE,
stars = TRUE, statistic = "std.error", # "[{conf.low} - {conf.high}]",
coef_omit = "SD|Cor",
coef_rename = c(
"poly(Crowdedness, 2)1" = "Crowdedness",
"poly(Crowdedness, 2)2" = "Crowdedness<sup> 2</sup>",
"Typologygreen" = "Green",
"Familiarity" = "Familiarity",
"poly(Crowdedness, 2)1Typologygreen" = "Crowdedness x Green",
"poly(Crowdedness, 2)2Typologygreen" = "Crowdedness<sup> 2</sup> x Green",
"poly(Structure, 2)1" = "Structure",
"poly(Structure, 2)2" = "Structure<sup> 2</sup>",
"poly(Structure, 2)1Typologygreen" = "Structure x Green",
"poly(Structure, 2)2Typologygreen" = "Structure<sup> 2</sup> x Green",
"poly(Scenic, 2)1" = "Scenic",
"poly(Scenic, 2)2" = "Scenic<sup> 2</sup>",
"poly(Scenic, 2)1Typologygreen" = "Scenic x Green",
"poly(Scenic, 2)2Typologygreen" = "Scenic<sup> 2</sup> x Green",
"poly(Width, 2)1Typologygreen" = "Width x Green",
"poly(Width, 2)2Typologygreen" = "Width<sup> 2</sup> x Green"
),
title = "Non-linear mixed effects models of Perceived Beauty \\label{tab:models_beauty}",
# notes = "PC (log) : Average pedestrian count per stimulus" ,
output = "output/table_models_perceived_beauty.tex"
)
3.4 Determinants of Psychological Reaction
rm(preds, dat, mods, vars)
preds <- data.frame()
dat <- data.frame()
mods <- list()
outcome_vars <- c("WillingnessToWalk", "Valence", "Arousal", "Restoration", "Presence")
for (i in 1:length(outcome_vars)) {
resp <- outcome_vars[i]
print(resp)
f <- paste(resp, " ~ Typology / (poly(Crowdedness, 2) * Beauty) + Familiarity + (1|Participant) + (1|Video)")
f
model <- glmmTMB::glmmTMB(as.formula(f), data = na.omit(df[c(resp, "Beauty", "Typology", "Crowdedness", "Familiarity", "Participant", "Video")]))
mods[[i]] <- model
param <- parameters::parameters(model, effects = "fixed")
# Predictions
preds <- estimate_relation(model, at = c("Typology", "Crowdedness", "Beauty"), length = c(50, 7)) |>
mutate(Outcome = resp) |>
rbind(preds)
# Data for density
dat <- rbind(
dat,
data.frame(
Typology = df$Typology,
Crowdedness = df$Crowdedness,
video_mean_pedcounts = df$video_mean_pedcounts,
Predicted = df[[resp]],
Outcome = resp
)
)
# Random (simplified model)
# f <- paste(resp, " ~ mean_pedcounts + Beauty + (mean_pedcounts + Beauty |Participant)")
# model <- glmmTMB::glmmTMB(as.formula(f), data = df)
# # dfsub <- modelbased::estimate_grouplevel(model) |>
# modelbased::reshape_grouplevel(indices = "Coefficient") |>
# select(c(1, 3, 4)) |>
# group_by(Participant) |>
# slice(1) |>
# data_addprefix(paste0(resp, "_"), exclude = "Participant") |>
# merge(dfsub, by = "Participant")
}
## [1] "WillingnessToWalk"
## [1] "Valence"
## [1] "Arousal"
## [1] "Restoration"
## [1] "Presence"
names(mods) <- c("Willingness to Walk", "Valence", "Arousal", "Restoration", "Presence")3.4.1 plots psych
dat$Outcome <- recode(dat$Outcome, "WillingnessToWalk"="Willing.-to-walk")
p_psych <- preds |>
dplyr::filter(Predicted <= 7) |>
mutate(Beauty = fct_rev(as.factor(Beauty)),
Outcome = recode_factor(Outcome, "WillingnessToWalk"="Willing.-to-walk"),
Outcome = fct_relevel(Outcome, c("Willing.-to-walk", "Valence", "Arousal", "Restoration", "Interest", "Presence"))) |>
ggplot(aes(x = Crowdedness, y = Predicted)) +
stat_density_2d(data = dat, aes(fill = ..density..), geom = "raster", contour = FALSE) +
geom_line(aes(size = Beauty, color = Typology, group = interaction(Beauty, Typology))) +
scale_size_manual(values = seq(from = .05, to = 2, length.out= 7), breaks= c(1:7)) +
scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) +
scale_alpha_ordinal(range = c(1, 0.1)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), breaks = c(1:7)) +
theme_modern(axis.title.space = 10) +
theme(
strip.placement = "outside",
strip.text.x = element_text(hjust = 0, size = 12),
axis.title.y = element_text(face = "bold"),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
labs(x = "Subjective Crowdedness", y = "Psychological Outcome", title = "Psychological effect of Crowdedness and its modulation by Beauty") +
ggnewscale::new_scale_fill() +
scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
ggside::geom_xsidedensity(data=dat, aes(y=stat(density), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
ggside::geom_ysidedensity(data=dat, aes(x=stat(density), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
ggside::theme_ggside_void() +
facet_wrap(~Outcome, strip.position = "left", scales = "free")
p_psych
3.4.2 Table output
modelsummary(mods, stars = TRUE)| Willingness to Walk | Valence | Arousal | Restoration | Presence | |
|---|---|---|---|---|---|
| (Intercept) | 0.078 | 1.675*** | 2.893*** | −0.320*** | 2.490*** |
| (0.053) | (0.044) | (0.067) | (0.064) | (0.070) | |
| Typologygreen | 0.064 | 0.028 | 0.007 | 0.540** | 0.062 |
| (0.155) | (0.127) | (0.196) | (0.173) | (0.138) | |
| Familiarity | 0.052*** | 0.066*** | 0.013 | 0.013+ | 0.151*** |
| (0.007) | (0.006) | (0.009) | (0.008) | (0.006) | |
| Typologyurban × poly(Crowdedness, 2)1 | −5.507 | −8.341** | 12.964** | 0.898 | 15.066*** |
| (3.460) | (2.839) | (4.398) | (3.908) | (3.062) | |
| Typologygreen × poly(Crowdedness, 2)1 | −17.504 | −36.489*** | −1.308 | −19.692 | 15.132 |
| (13.306) | (10.936) | (16.952) | (15.026) | (12.081) | |
| Typologyurban × poly(Crowdedness, 2)2 | −26.159*** | −29.599*** | 12.199** | −16.120*** | 3.829 |
| (3.361) | (2.753) | (4.265) | (3.777) | (3.071) | |
| Typologygreen × poly(Crowdedness, 2)2 | −25.949* | −40.108*** | −3.835 | −22.872+ | −16.571 |
| (12.189) | (9.979) | (15.459) | (13.692) | (11.337) | |
| Typologyurban × Beauty | 0.758*** | 0.638*** | 0.206*** | 0.799*** | 0.255*** |
| (0.008) | (0.007) | (0.010) | (0.010) | (0.007) | |
| Typologygreen × Beauty | 0.741*** | 0.636*** | 0.083* | 0.763*** | 0.250*** |
| (0.026) | (0.021) | (0.032) | (0.028) | (0.023) | |
| Typologyurban × poly(Crowdedness, 2)1 × Beauty | 0.625 | 0.952 | 6.878*** | −1.104 | −2.505*** |
| (0.731) | (0.599) | (0.930) | (0.824) | (0.662) | |
| Typologygreen × poly(Crowdedness, 2)1 × Beauty | 0.292 | 4.187* | 4.824+ | 1.049 | −2.959 |
| (2.267) | (1.861) | (2.883) | (2.552) | (2.087) | |
| Typologyurban × poly(Crowdedness, 2)2 × Beauty | 3.214*** | 4.649*** | −3.547*** | 1.357+ | −0.503 |
| (0.723) | (0.591) | (0.917) | (0.811) | (0.667) | |
| Typologygreen × poly(Crowdedness, 2)2 × Beauty | 2.136 | 5.612*** | −0.186 | 2.660 | 1.931 |
| (2.072) | (1.696) | (2.629) | (2.327) | (1.931) | |
| SD (Intercept Participant) | 0.528 | 0.461 | 0.670 | 0.742 | 1.107 |
| SD (Intercept Video) | 0.175 | 0.151 | 0.229 | 0.215 | 0.000 |
| SD (Observations) | 1.010 | 0.825 | 1.280 | 1.129 | 0.951 |
| Num.Obs. | 10899 | 11009 | 11009 | 11107 | 11106 |
| R2 Marg. | 0.537 | 0.554 | 0.119 | 0.486 | 0.239 |
| R2 Cond. | 0.645 | 0.669 | 0.326 | 0.650 | |
| AIC | 32248.3 | 28168.2 | 37791.7 | 35523.6 | 31831.3 |
| BIC | 32365.0 | 28285.1 | 37908.6 | 35640.7 | 31948.3 |
| ICC | 0.2 | 0.3 | 0.2 | 0.3 | |
| RMSE | 0.99 | 0.80 | 1.25 | 1.10 | 0.94 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
modelsummary(mods,
stars = TRUE, statistic = "[{conf.low} - {conf.high}]", escape = T,
coef_omit = "SD|Cor",
coef_rename = c("Typologygreen"="Green",
"Familiarity" = "Familiarity",
"Typologyurbanpoly(Crowdedness, 2)1"="Urban * Crowdedness",
"Typologygreenpoly(Crowdedness, 2)1"="Green * Crowdedness",
"Typologyurbanpoly(Crowdedness, 2)2"="Urban * Crowdedness<sup> 2</sup>",
"Typologygreenpoly(Crowdedness, 2)2"="Green * Crowdedness<sup> 2</sup>",
"TypologyurbanBeauty" = "Urban * Beauty",
"TypologygreenBeauty" = "Green * Beauty",
"Typologyurbanpoly(Crowdedness, 2)1Beauty" = "Urban * Crowdedness * Beauty",
"Typologygreenpoly(Crowdedness, 2)1Beauty" = "Green * Crowdedness * Beauty",
"Typologyurbanpoly(Crowdedness, 2)2Beauty" = "Urban * Crowdedness<sup> 2</sup> * Beauty",
"Typologygreenpoly(Crowdedness, 2)2Beauty" = "Green * Crowdedness<sup> 2</sup> * Beauty"
),
title = "Non-linear mixed effects models of psychological reaction\\label{tab:models_reaction}",
notes = "PC (log) : Average pedestrian count per stimulus" ,
output = "output/table_models_psych_reaction.tex"
)3.5 joint plots
p_reg <- ((p_crowd + guides(fill = "none", color = "none") | p_beauty + theme(legend.position = "none")) / p_psych + theme(legend.position = 'bottom') ) +plot_annotation(tag_levels = "A")
p_reg
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 420 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 15 rows containing non-finite values (`stat_density()`).
## Warning: Removed 405 rows containing non-finite values (`stat_density()`).
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
ggsave(p_reg, filename = "output/plots_reg.jpg", width = 13, height = 13, dpi = 300)
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 420 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 15 rows containing non-finite values (`stat_density()`).
## Warning: Removed 405 rows containing non-finite values (`stat_density()`).
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.pvalue <- function(x, ...) {
# Construct vectors of data y, and groups (strata) g
y <- unlist(x)
g <- factor(rep(1:length(x), times=sapply(x, length)))
if (is.numeric(y)) {
# For numeric variables, perform a standard 2-sample t-test
# p <- t.test(y ~ g)$p.value
p <- BayesFactor::ttestBF(x=y, y= g)
p <- p@bayesFactor$bf
# print(p)
} else {
# For categorical variables, perform a chi-squared test of independence
p <- chisq.test(table(y, g))$p.value
}
# Format the p-value, using an HTML entity for the less-than sign.
# The initial empty string places the output on the line below the variable label.
c("", sub("<", "<", format.pval(p, digits=3, eps=0.001)))
}
library(kableExtra)
t1 <- table1::table1(~ Valence + Arousal + Restoration + Beauty + `Willingness to walk` + Interest + Familiarity + Structure + Presence + Crowdedness | Typology, overall = F,
data=dfenv, extra.col=list(`BF`=pvalue))
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
t1| Urban (N=408) |
Green (N=72) |
BF | |
|---|---|---|---|
| Valence | |||
| Mean (SD) | 4.75 (0.638) | 5.81 (0.525) | 1175 |
| Median [Min, Max] | 4.81 [2.53, 6.19] | 5.94 [3.36, 6.71] | |
| Arousal | |||
| Mean (SD) | 3.90 (0.547) | 3.29 (0.546) | 1018 |
| Median [Min, Max] | 3.91 [2.14, 5.33] | 3.33 [1.97, 4.46] | |
| Restoration | |||
| Mean (SD) | 3.25 (0.817) | 4.85 (0.532) | 606 |
| Median [Min, Max] | 3.23 [0.870, 5.36] | 4.95 [2.91, 5.74] | |
| Beauty | |||
| Mean (SD) | 4.41 (0.761) | 5.89 (0.514) | 958 |
| Median [Min, Max] | 4.36 [1.96, 6.41] | 6.02 [4.27, 6.83] | |
| Willingness to walk | |||
| Mean (SD) | 3.63 (0.750) | 4.82 (0.584) | 787 |
| Median [Min, Max] | 3.68 [1.30, 5.44] | 4.98 [2.52, 5.74] | |
| Interest | |||
| Mean (SD) | 4.33 (0.832) | 5.35 (0.608) | 939 |
| Median [Min, Max] | 4.33 [2.22, 6.24] | 5.40 [3.80, 6.23] | |
| Familiarity | |||
| Mean (SD) | 4.10 (0.559) | 4.49 (0.602) | 1130 |
| Median [Min, Max] | 4.13 [2.05, 5.68] | 4.52 [2.82, 5.76] | |
| Structure | |||
| Mean (SD) | 4.37 (0.989) | 5.48 (0.606) | 847 |
| Median [Min, Max] | 4.45 [1.40, 6.54] | 5.61 [3.29, 6.32] | |
| Presence | |||
| Mean (SD) | 4.23 (0.332) | 4.71 (0.339) | 1419 |
| Median [Min, Max] | 4.25 [2.91, 5.13] | 4.77 [3.82, 5.48] | |
| Crowdedness | |||
| Mean (SD) | 4.42 (1.42) | 2.99 (1.47) | 509 |
| Median [Min, Max] | 4.56 [1.22, 6.88] | 2.84 [1.14, 5.96] |
library(ggdist)
penv <- dfenv %>%
pivot_longer(Restoration : Crowdedness) %>%
ggplot(aes(x = name, y = value, colour = Typology)) +
# geom_boxplot( position = position_dodge()) +
# geom_line(aes(group = name_df_keep), alpha =.2)+
geom_jitter(position = position_jitterdodge(), alpha = 0.5) +
ggdist::stat_halfeye(position = position_dodge(), aes(fill = Typology),
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
## move geom to the right
justification = -.4,
## remove slab interval
.width = 0,
point_colour = NA) +
scale_color_manual(values = c("Urban" = col_urban, "Green" = col_green)) +
scale_fill_manual(values = c("Urban" = col_urban, "Green" = col_green)) +
theme_classic() +
# gghalves
labs(y = "Rating", x = "Dimension")
penv